
#ifndef __IO_OFFS_H__
#define __IO_OFFS_H__

#include <string>
#include <fstream>
#include <vector>
#include <cassert>
#include <set>
#include <cmath>

using namespace std;

extern vector<Point3>      mesh_pos;
extern vector<Vector3>     mesh_normals;
extern vector<unsigned>    mesh_faces;
extern int                 nfaces;
extern int                 nvertices;

/// Read the OFF model
/// @arg filename the file name
/// @return true if it succeed
bool readOFFfile(const char *filename) {
   
   ifstream filein(filename);
   
   if (not filein.is_open()) {
      cout<<"I/O error: the file doesn't exist "<<endl;
      return false;
   }
   else {
      string type;
      int nt; float x, y, z; unsigned nv, v0,v1,v2;
      filein>>type;
      filein>>nvertices>>nfaces>>nt;
      assert(type == "OFF");
      assert(nvertices>0 and nfaces>0);
      
      mesh_pos.resize(nvertices);
      mesh_normals.resize(nvertices);
      mesh_faces.resize(nfaces*3);
      
      /// load the mesh positions
      for (int i=0; i<nvertices; ++i) {
         filein>>x>>y>>z;
         mesh_pos[i] = Point3(x,y,z);
      }
      
      /// load the array of faces
      std::vector<std::set<int> > facespervertex(nvertices);
      for (int i=0; i<nfaces; ++i) {
         filein>>nv>>v0>>v1>>v2;
         mesh_faces[i*3+0] = v0;
         mesh_faces[i*3+1] = v1;
         mesh_faces[i*3+2] = v2;
         facespervertex[v0].insert(i);
         facespervertex[v1].insert(i);
         facespervertex[v2].insert(i);
      }
      
      /// compute one normal per face
      std::vector<float> facenormals(3*nfaces);
      for (int i=0; i<nfaces; ++i) {
         const Point3 &p0 = mesh_pos[mesh_faces[i*3+0]];
         const Point3 &p1 = mesh_pos[mesh_faces[i*3+1]];
         const Point3 &p2 = mesh_pos[mesh_faces[i*3+2]];
         float A[3] = {p1[0]-p0[0], p1[1]-p0[1], p1[2]-p0[2]};
         float B[3] = {p2[0]-p0[0], p2[1]-p0[1], p2[2]-p0[2]};
         
         /// Computes the cross product: A x B
         facenormals[i*3+0] = A[1]*B[2] - A[2]*B[1];
         facenormals[i*3+1] = A[2]*B[0] - A[0]*B[2];
         facenormals[i*3+2] = A[0]*B[1] - A[1]*B[0];
      }
      
      /// compute one normal per vertice
      for (unsigned int i=0; i<facespervertex.size(); ++i) {
         const std::set<int> &incidentfaces = facespervertex[i];
         std::set<int>::const_iterator it = incidentfaces.begin();
         float s[]={0.0,0.0,0.0}; int c=0;
         while (it != incidentfaces.end()) {
            const int &f = *it++; ++c;
            s[0] += facenormals[f*3+0];
            s[1] += facenormals[f*3+1];
            s[2] += facenormals[f*3+2];
         }
         s[0] /= (float)c;
         s[1] /= (float)c;
         s[2] /= (float)c;
         float sl = sqrt(s[0]*s[0] + s[1]*s[1] + s[2]*s[2]);
         if (sl > 0) mesh_normals[i] = Vector3(s[0],s[1],s[2])/sl;
         else        assert(false);
      }
      
      return true;
   }
}

/// Save the OFF model
/// @arg filename the file name
/// @return true if it succeed
bool saveOFFfile(const char *filename) {
   
   ofstream fileout(filename);
   if (not fileout.is_open()) {
      cout<<"I/O error "<<endl;
      return false;
   }
   else {
      assert(not mesh_pos.empty());
      assert(not mesh_faces.empty());
      assert(nfaces>0 and nvertices>0);
      assert(mesh_pos.size() == (unsigned)nvertices and mesh_faces.size()/3 == (unsigned)nfaces);
      
      fileout.precision(10);
      fileout.setf(ios::fixed, ios::floatfield);
      fileout<<"OFF"<<endl;
      fileout<<nvertices<<" "<<nfaces<<" "<<"0"<<endl;
      for (unsigned i=0; i<mesh_pos.size(); ++i) 
         fileout<<mesh_pos[i]<<endl;
      for (unsigned i=0; i<mesh_faces.size(); i+=3) 
         fileout<<"3 "<<mesh_faces[i+0]<<" "<<mesh_faces[i+1]<<" "<<mesh_faces[i+2]<<endl;
      return true;
   }
}

#endif
